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ANALYSIS BY CONFORMAL MAPPING OF STEADY FROZEN-LAYER 
CONFIGURATION INSIDE COLD RECTANGULAR CHANNELS 
CONTAINING WARM FLOWING LIQUIDS 
by Robert Siegel and Joseph M. Savino 
Lewis Research Center 

SUMMARY 

Conformal mapping was applied to determine the shape of the solidified region that 
forms inside a cooled rectangular channel containing flowing liquid. The channel walls 
are at a uniform temperature below the freezing temperature of the liquid. The liquid 
supplies energy by convection to the frozen region interface, and the shape of the inter- 
face adjusts so that this energy can be conducted through the region to the cold walls. 

The interface shape is obtained by mapping the solidified region into a rectangle in the 
potential plane and then determining the conformal transformations between the potential 
and physical planes . For a given aspect ratio of the duct, the configuration of the frozen 
region and the heat flow through it depend on a single physical parameter involving the 
liquid heat-transfer coefficient, wall and liquid temperatures, and the freezing tempera- 
ture. For a given value of this physical parameter in a certain range, it was mathemat- 
ically found that there can be two shapes of the frozen region, although only one of these 
shapes would be physically stable. When the cooling parameter is above a certain value 
for each aspect ratio, the duct will freeze completely. 


INTRODUCTION 

The aim of this work is to study the two-dimensional solidification of warm liquids 
that flow through cooling channels of heat exchangers. Some compact heat exchangers 
utilize straight flow passages of a rectangular cross section, and in certain applications 
when warm liquids are being cooled, the conditions of temperature and flow can cause 
the channel walls to cool below the freezing temperature of the liquid. Consequently, 
some of the liquid freezes to form a frozen region between the cold walls and the flowing 
liquid. This provides an insulating region that reduces the amount of energy that can be 
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transferred from the warm liquid to the walls. The heat flow through the frozen region 
must be computed to evaluate the system performance. 

Little work has been done previously on freezing inside two-dimensional geometries. 
In reference 1, transient freezing inside a square prism was analyzed. The liquid was 
at the fusion temperature so that no convective energy was being supplied to the frozen 
interface. However, for the transient situation, energy is being released at the interface 
as a result of the heat of fusion, and this energy must be conducted away through the 
frozen region. Several references (refs. 2 to 5) present analyses of solidification in cir- 
cular tubes or the gap between parallel plates, which are treated as one- dimensional 
situations at each axial location. The behavior in a circular tube is used later in the 
section RESULTS AND DISCUSSION to aid in the discussion of the present results. 

The heat conduction solution in a two-dimensional region with known boundaries 
would ordinarily not present a great difficulty because numerical solutions of Laplace's 
equation can be utilized if an analytical solution is not available. In the present situa- 
tion, however, the shape of the frozen region is unknown. The shape adjusts so that the 
energy supplied by convection from the flowing liquid to the frozen interface can pass by 
conduction through the solid to the cold walls. Thus, the frozen shape must be found to 
compute the heat being transferred. 

The shape of the frozen region and the heat transfer through it is found by applying a 
conformal mapping technique that was developed in reference 6. The conformal mapping 
technique is feasible because the boundaries of the frozen region are at constant temper- 
atures. The liquid - frozen- layer interface is at the equilibrium freezing point while the 
channel walls are at a specified constant temperature below the freezing point. As a re- 
sult, the frozen region can be represented as a rectangle in the potential plane. The 
rectangle can be mapped conformally into the physical plane to determine the frozen con- 
figuration. 

In a heat exchanger passage, the liquid temperature would be changing as it flows 
along the passage length. The present analysis is two-dimensional and, hence, does not 
include effects such as axial conduction. Usually, in a channel, the heat conduction in 
the axial direction is small compared with that transverse to the fluid flow direction. 
Hence, the present results can be applied locally along the channel length to compute the 
frozen- layer development in a channel. 


SYMBOLS 


A dimensionless half -width, a /y 
a half-width of long side of channel 
B dimensionless half -width, b/y 
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b 


half-width of short side of channel 


Cf, Cg, Cg 
c,d 

F 

h 


constants 

intermediate parameters in mapping 


elliptic integral of the first kind, F (<p, k) = 


f 9 d <P 

Vl - k 2 sin 2 <p 


heat-transfer coefficient 


K 

k 

N 

n 

Q 

q 

r 

s 

T 

t 

u 

W 

X,Y 

x,y 

z 

y 

C 




complete elliptic integral of the first kind, K(k) = 



d<p 


V 1 - k 2 sin 2 <p 


thermal conductivity of solidified material 
dimensionless outward normal distance, n/y 
outward normal distance from interface 

heat flow rate through frozen region per unit length of channel 

heat flow per unit area and time 

radius 

length of frozen interface 
dimensionless temperature, (t - t w )/(t^ - t w ) 
temperature; dummy integration variable 
mapping variable, £ + it] 
analytic function, -T + ii// 
dimensionless coordinates, x/y, y/y 
Cartesian coordinates in physical plane 
dimensionless physical plane, X + iY 


length scale parameter, 

quantity defined as - — 

9X 


k % ~ *w 
t l " *f 

+ i ^ 

9Y 


coordinates of u plane 
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6 argument of £ 

x{/ imaginary part of W 

o> complex quantity, In | £ | + id 

Subscripts: 

f at freezing point 

i inner 

l liquid 

o outer 

s at frozen interface 
w wall 


ANALYSIS 

Solution for General Case 

The model of a typical cross section of a rectangular flow channel in which a steady- 

state frozen layer has formed is shown in figure 1. The warm liquid is at a constant 

bulk temperature t^, and the local heat-transfer coefficient around the interface h^ is 

assumed constant. The channel walls are assumed to be held at a constant temperature 

t w , which is below the liquid freezing temperature t f . The warm liquid is supplying the 

heat flux h^(t^ - t f ) to the liquid - frozen- layer interface. The coordinates x g and y g 

of the interface are unknown and are dependent on the duct aspect ratio a/h, the freezing 

temperature t f , the thermal conductivity k of the frozen material, and the parameters 

hj, tj, and t w . The objectives of this analysis are to determine the coordinates x g 

and y of the interface and the heat transfer Q to the channel walls per unit channel 
•'s 

length. 

Basically, the solution to the problem involves solving Laplace's heat conduction 
equation in the frozen region subject to the known boundary conditions. That is, in the 
frozen region 

A + ift = 0 (l) 

3x 2 3y 2 


4 



On the interface (x 0 , y_), 
s s 


kii 

Sn 


- t f ) 


x s>y s 


t (x s , y s ) = t f 


( 2 ) 


On the channel walls, 


t = = constant 


(3) 


Because of symmetry, only one quadrant of the channel cross section need be considered 
(fig. 2(a)). The symmetry provides the two conditions 


at 

0X 

at 

3y 


= o 


x=a 


= 0 


y=b 


(4) 


The equations are now normalized by letting 


T s. 


t - t 


w 





where 



r 55 


*<*1 - v> 

h i^i ~ v 


The normalized Laplace equation and boundary conditions take the form 
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(5) 




= 0 


On the interface (X„, Y_) (fig. 2(b)), 
s s 


3T 

3N 


= 1 


X s’ Y s 


T ( X s’V =1 ^ 


( 6 ) 


On the walls at X = 0, 0 < Y < B and at 0 < X < A, Y = 0, 


T = 0 


(7) 


On the symmetry lines, 


3T 

ax 

3T 

aY 




= 0 


X=A 


= 0 


Y=B 


(8) 


The conformal mapping method developed in reference 6 is applied herein. Let 


W = -T + ixp (9) - 

where W is the complex potential having -T and \p as the real and imaginary parts. 

The functions -T and xf/ are conjugate harmonic functions of X and Y, and, there- 
fore, W is an analytic function of the complex variable Z. The coordinates (X„, Y_) of 
the interface and the total heat transfer Q are found by representing the frozen region 
in the W plane and then finding a conformal transformation between the W and Z 
planes. 

To find the relation between quantities in the W and Z planes, note that the de- 
rivative of an analytic function is independent of direction. Hence, 

^=-9 1 + 1 M (1 0) 

dz ax ax 
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The Cauchy-Riemann equations can be used to relate the derivatives of the real and 
imaginary parts of W, giving 


3i p _ 3T 
3X 3Y 


Using equation (11) to eliminate dip/dX from equation (10) gives 


dW _ aT [ ■ 3T 
dZ 3X 3Y 


( 11 ) 


(12) 


Now define £ as 


?s_iZ + ilZ ( 13 ) 

ax 3Y 

so that equation (12) becomes 

dW _ ^ 
dZ 

Integrating gives 

Z = f - dW + constant (14) 

J C 

Equation (14) shows that the connection between Z and W involves an integration 
containing £; hence, £ must be related to W before the integration can be performed. 
The general shapes of frozen region in the W and £ planes are known from the 
physical conditions of the problem, as shown in the next section, and the relation be- 
tween them is found by conformal mapping. 

Conformal mapping transformations . - Refer now to figures 2(b) and 3. Since the 

wall 6712^ and the interface 345 are isotherms, they become vertical lines in figure 3(a). 
No heat flows across boundaries 23 and 56; hence, these boundaries are along the direc- 
tion of heat flow. The heat flow lines are normal to the isotherms, so these boundaries 
form the top and bottom of a rectangle in figure 3(a). 

The frozen region in the temperature derivative plane is shown in figure 3(b). Be- 
cause 3T/3N = 1 on 345, this boundary becomes a part of the unit circle in figure 3(b). 
Along 123, 3T/3X = 0, while along 567, 3T/3Y = 0 (fig. 2(b)), so that the frozen region 
maps into a quarter of the unit circle in figure 3(b). If a transformation between the re- 
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gions in figures 3(a) and (b) can be found, the integration in equation (14) can be per- 
formed, and the line 345 in figure 3(a) can be transformed into the physical plane to ob- 
tain the shape of the frozen interface. Several successive intermediate transformations 
are necessary to relate £ and W. The derivation of the necessary transformations is 
now presented. 

The region in the derivative plane £ (fig. 3(b)) is mapped into the region of the co 
plane shown in figure 3(c) by the following transformation: 

w = In £ = ln|£| +i0 (15) 

The generalized rectangular region of the co plane is mapped into the upper half of 
the u-plane in the manner indicated by figures 3(c) and (d) by applying the Schwarz- 
Christoffel transformation 


dw C 1 



(16) 


Integrating equation (16) gives 


w = Cj 


In 



+ C 


2 


The constants C ^ = -1/2 and C 2 = iff are chosen so that the points 3 and 5 of the u 
plane map into the points 3 and 5 of the co plane such that the horizontal lines are 
located in the w plane at 8 = n and n/2. Therefore, 



(17) 


Now the transformations relating ( to W are completed by mapping the upper half 
of the u plane into the rectangle in the W plane by applying the Schwarz-Christoffel 
transformation in the following form: 


dW ^3 

Vc u 2 - l)(u + c)(u - d) 


(18) 


The boundary of the rectangle in the W plane maps into the horizontal axis of the u 
plane shown in figure 3(d). 

Integration into physical plane. - With the equations (15), (17), and (18) available, 
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the integration in equation (14) can be carried out. The integration is performed through 
the use of the intermediate variable u 


Z = 


fi 


dW 
£(u) du 


du + constant 


(19) 


To determine £(u), solve equation (15) for £ to obtain £ = e w and substitute c o from 
equation (17) to obtain (note that e 17r = -1) 


?(u) 




Substitute equations (18) and (20) into equation (19) to obtain 


( 20 ) 


z = -c. 


/k 


u + Vu 2 - 1 


- 1 1/2 


l)(u - l)(u + c)(u - d) 


du + constant 


( 21 ) 


The numerator of equation (21) can be expressed in the form 




as can be verified by squaring both sides. Then equation (21) becomes 


Z = 


Vi 


^ du + f du_ 

J V(u - l)(u + c)(u - d) J V(u + l)(u + 


c)(u - d) 


+ constant (22) 


On the frozen interface, u = £ + iO for -1 < | < 1. Hence, along the interface, 


Z c - Z- = 
S 5 


Vi 


/' 

•a 


d| 


V(i - m + c)(d - ?) 


+ i f 4S- 

J x Vd + m + 


c)(d - |) 


for -1 < | < 1 (23) 
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To find Cg, use the condition Z g - Z ? = iB in equation (22) evaluated along the 
real axis of the u-plane from £ = d to £ = °° to obtain 



The integrals in equation (24) are complete elliptic integrals and can be written as 



By using equation (25), the Cg is eliminated from equation (23) to give 



The Zg will be eliminated by first applying equation (26) at point 3, which is at 
£ = -1, to give 



Now, from figure 2(b), Zg = A + iB - Re (Zg - Z g ). By use of equation (27) to obtain 
Re (Zg - Zg), Zg is eliminated from equation (26) to give 



■ ii 





To obtain X , take the real part of equation (28),, which gives 


X s _ A _ Vc + d 
B B 


r — a. 

75 J. i Vd + m + 


c)(d - 4) 




for -1 < | < 1 


(29) 


K 


Similarly, Y_ is found by taking the imaginary part of equation (28) 

S 


/ 

1 _ V c + d 


d|_ 


Vd - m + c)(d - 4) 


▼ c + dy \Jc + d 


for -1 < 4 < 1 


(30) 


K 


The integrals in the numerators of equations (29) and (30) can be transformed into 
elliptic integrals of the first kind to give 


B 


A 

B 


in“l 7(c + d)(l + 4) t /l + 

V( l + d)(c+ 4)’ V C + 

\l c + dy c + d, 


for -1 < 4 < 1 


(31) 
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3!.. 

B 


(sin' 1 JHWH, JIH' 

[_ f (1 + c)(d - |) T c + d_ 


K 


(\S) +K (^ 


for -1 < ^ < 1 


(32) 


The aspect ratio A/B must be expressed in terms of the quantities involved in the 
mapping. The length A = Zg - Z^, where Zg and Zj correspond to £ = -c and -«>. 

Then, from equation (21) (note that in this interval of £, )[v? - 1 = ^/e i7r | £+ 1 |e iir | £ - 1 1 

= e ilT ^~- 1 = -Vl 2 - D. 


A = -C, 


/> 


Vi 2 - 1 


D(« - D(S + c)(| - d) 


1/2 




Eliminate Cg by using equation (25) to obtain 


A=-i 

B 


Vc 4- d 

Vi 


/T- 
— / (« + 
+ d t/ - 00 L- 




1)(« - D(« + c)(« - d) 


1/2 




K 


(VS) ^ k (vS ’ 


Change the variable in the numerator by letting £ = -t to yield, 


_ a _ _ V< 


A_ a 
B b 


: + d 




(t+l)(t- l)(t- c)(t + d) 


Vi 


K 


2 

* JjnYJjiii 

\jld + cJ yTd+cy 

( 0 )+ k ( 0 ) k (VS) +k (v§J 


(33) 


Equation (33) relates the aspect ratio of the duct to the quantities c and d that are 
involved in the mapping. 

Heat flow through frozen regio n. - For completion of the solution, the heat flow 
through the frozen region needs to be determined. Also, the mapping parameters 
c and d need to be related to the basic physical variables in the problem, such as the 
convective heat-transfer coefficient, the liquid temperature, and the freezing tempera- 
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ture. The heat flow rate per unit channel length through the frozen region and into the 
cold wall is four times the value for one quadrant; that is, 


Z y=b /»x=a 

k ! dy+ [ 

_ • x_o «4o 


kit 

‘by 


dx 


y=o 


Normalizing gives 


4k(t. 


Apply the Cauchy-Riemann conditions 


/ Y=B /»X=A 

il dY + / il 

.J 3 X x=o 4 =0 9Y 


dY 


|Y=0 


_ 9T _ 9^ . 0T _ jh£ 

ax 8Y ’ aY ax 


to obtain 


Q 


4k(t f 


/ Y=B 
_ J 


3Y 


dY + 


X=0 


✓»X=A 

*4=o 


dip 

ax 


dX 


Y=0 


or 


Q 


= ~ip(0, B) + ip( A,0) = (-i)[w(2) - W(6)] = (-i)[w(3) - W(5)] 


4k(t f - t w ) 

The foregoing can be evaluated as 

Q 


4k (tf - t w ) 

which can also be written as 

Q_ 

4k(t f - 


= -1 


/ 


point 3 


dW 

du 


du 


point 5 


V-/; 


dW 

du 


du 
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Substitute dW/du from equation (18) and note that u = £ on the real axis 


q = iC r ^ 

4k ^f ■ V *4=-l V(4 2 - m + c )(| - d) 

The constant Cg is conveniently foimd here by noting from figure 3(a) that 
1 = W(6) - W(5), which can be written as 


( 34 ) 


/ d /»d 

™ du.C, / d * — 

du Vd 2 - Du + cMi - d) 


Solve for C 


3 


C 3 " 


1 f dg 

*4=i V(^ 2 - i)u +c)(d- k) 


(35) 


Substitute equation (35) in (34) to obtain 


Q 


i: 






na + c)(d - o 


4k(t f - t w ) />d 


/' 


d* 




(r - i)(s + c)(d - 1) 


(36) 


The integrals are complete elliptic integrals, and equation (36) can be put in the form 


Q 

4k(t f - t^) 


K 


K 


J 2(c + 

_f (c + 1)(< 

m 


d) 


l)(d + 1)_ 

- l)(d - 1) 
l)(d + 1) 


(37) 
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Relation between phys ical conditions and mapping . - Finally, an expression is 
needed to relate the parameter 


— = k *f ~ *w 

B V h ~ h 

which contains the controllable variables, to the mapping quantities c and d. This 
relation can be found by noting that two expressions have been obtained for Cg, equa- 
tions (25) and (35). Equating these two expressions gives 



The numerator can be expressed as a standard elliptic integral giving 


k *f ~ V 
h £ b *7 " *f 


/ 2(c + d) k 

!(c - l)(d - 1) 

V(c + l)(d+ 1) |V 

(c + l)(d + 1)_ 



(38) 


Evaluation of solution . - The solution to the problem consists of equations (31), (32), 
(33), (37), and (38). In using this solution, the aim is to calculate the total heat trans- 
fer Q and the coordinates x_,y e of the liquid - frozen-layer interface, given a known 
set of values for t w , t^, t^, h^, k, a, and b. 

The first step is to evaluate equation (33) for a variety of c's and d's to find a set 
of pairs of c and d that results in the specified aspect ratio a/b. Then, each pair of 
c and d is used in equations (37) and (38) to determine the values of Q/4k(t^ - t w ) and 
(k/h^b) [(tf - t w )/(t j - t^)j. The c and d pair and the corresponding aspect ratio are 
used in equations (31) and (32) to obtain the coordinates of the liquid - frozen-layer inter- 
face. By use of the various c and d pairs for a fixed aspect ratio, a family of inter- 
face shapes can be plotted, and the values of (k/h^b) [~(t f - t w )/(t^ - t f )] and Q/4k(t f - t w ) 
are known for each shape. The results of these calculations are given in the section 
RESULTS AND DISCUSSION. 
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Limiting Cases 


The previous analysis treated the general situation of a duct of any aspect ratio. 

Two limiting cases that are of interest and lead to simplified analytical results are 
shown in figure 4. 

In figure 4(a), the frozen region is sufficiently thin so that a one- dimensional region 
is formed in the central region along the long side. The following sections will show 
that this condition corresponds to having the mapping quantities c = 1 and d > 1. 

In figure 4(b), the frozen region is so thin that a one- dimensional region is formed 
in the central portion of each side. In this instance, the analysis will show that the 
mapping quantities have gone to the limits c = 1, d = 1. 

Limiting case for c - 1, d > 1: layer one- dime nsional in central region of long 
side . - With reference to figure 4(a), the frozen region becomes one- dimensional at 
the center of the long side. Increasing the aspect ratio will not change the shape of 
the frozen region but will only increase the extent of the one- dimensional region. 

Hence, a duct of infinite aspect ratio can be analyzed as shown in figure 5(a). 

The geometry is nondimensionalized in the same manner as for figure 2(b), which 
yields figure 5(b). The only difference in the two situations is that now points 2 and 3 
are both at infinity in the physical plane and, hence, also in the potential plane in fig- 
ure 6(a). The subsequent mappings as shown in figures 6(b), (c), and (d) then corre- 
spond to those in figures 3(b), (c), and (d). 

With points 2 and 3 at the same location in figure 6(d), the parameter c becomes 1. 
Equation (26) for coordinates on the frozen interface then becomes (note that K(0) = n/2) 


z s^- z s _vr^ 

B 2 



-1 < I < 1 (39) 


The distance Zg - Zg will be found to relate the interface to the coordinates of the duct. 
Using equation (22) with c = 1 and noting that Ug and u g correspond to | = 1 and 
4 = d in figure 6(d) result in 




1 

- m + D(d - T) 


i 

u + dVJTI 




1 < £ < d (40) 


16 



From equation (25), 



Eliminate Cg/i from equation (40) by using equation (41) and combine the result with 
equation (39) to eliminate Zg. These operations give 



Separate the real and imaginary parts to obtain 



The integrals in the numerator of X_/B can either be carried out or expressed as an 

s 

elliptic integral, and the integral in the numerator of the Y„/B expression can be 

s 
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expressed as an elliptic integral. These operations yield the expressions 




VdTl - Vd-l 
\Vd7i - VdTJ, 


+ln 

VdTT - Vd-T 




-1 < 5 < 1 (42) 



-1 < £ < 1 


(43) 


Equation (36) can be used to compute the heat flow by letting c = 1. Also, the inte- 
gral in the numerator is evaluated from £„ to 1, where | is the value of £ in 

a, a, 

equation (42) when X /B = A/B. This cuts off the infinite aspect ratio configuration 

s 

in figure 5(b) so that it extends only to the size of the finite aspect ratio duct. The 
heat flow is then 


Q 

4k (t f - V 



dj 

(J + 1)V(1 - g)(d - g) 

d| 

u + 1 )Vu - m - 1) 


Both integrals can be evaluated analytically, and after simplification, the heat flow be- 
comes 


Q 


= - In 


4k(t f - t w ) n 


2 V2(l + d) (d+ !)4 a + d + 4(d + 1) 


(1 + ^)(d - 1) 


d + 3 
d - 1 


( 44 ) 


Using equation (38) with c = 1 gives 
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1 


(45) 


k_ *f ~ *w _ 1 _ 



Frozen interface profiles are evaluated by choosing a value of d. The physical pa- 
rameter 1/B is found from equation (45). The corresponding heat flow and coordinates 
of the frozen region are found from equations (44), (42), and (43). 

Limiting case where c = 1 and d = 1: fro zen layer in corner region . - For a very 
thin frozen layer, the thickness is constant in the central portion of both sides of the duct 
away from the corners, as shown in figure 4(b). The configuration can then be found by 
joining four of the corner regions shown in figure 7(a). 

As in the preceding analyses, the quantities are nondimensionalized to yield figure 
7(b). The geometry is symmetric about the line 14. Since the pairs of points 2, 3 and 
5, 6 are at infinity, the mapping in the potential plane consists of two constant tempera- 
ture lines extending to infinity. The mappings, as shown in figure 8, are similar to 
those in figure 3. In the u plane, the pairs 5, 6 and 2, 3 will each be at a single point 
that can be fixed at (1, 0) and (-1, 0). Thus, the corner solution is a limiting case of the 
general solution for the situation when c = 1 and d= 1. Use equation (45) to eliminate B 
from equation (39) and then let d = 1 to obtain (note t is a dummy variable of integration) 


f r — h=- 1 — v= 

^ Jy l_(l+t)Vl-t (l-t)Vl+t_ 


dt 


At Z 4 , 4 = 0, so 


Z 4 - Z 5 = 




Then 


Z s- Z 4 = 


f '\— ! ._L_U 

J 0 ij 1 + t)Vn~t (i - 1) VTTt_ 


v5 /'[_] 

* J Q L(i + t)Vrri (i - t)Vm_ 


dt 


(46) 


Equation (46) is integrated to yield 
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z - Z. = 
s 4 


Vi - Vi - 1 Vi - Vi + I 

. i ( , n vs ♦ vr i , , v; ♦ vm 

Vi - 1 Vi- 1 


vi + 


Vi + 


( 47 ) 


The real part of Z g - Z 4 must approach Xg - X 4 as £ — 1, which gives 


X 5 -X 4 =Iln^-l 


77 VVi + 1/ 


(48) 


At a location far from the corner, the interface is flat and parallel to the wall. A heat 
balance yields 


«t f - v 


X c 


" " V 


so that 



Substituting into equation (48) and noting the symmetry of figure 7(b) gives 


X 4 = Y 4 = 1 + ± In 


*/i- 


(49) 


Now, the final expressions for the coordinates of the interface can be written. 

From symmetry, only the portion of the profile between Z 4 and Zg need be considered. 
Using equation (49) and the real and imaginary parts of equation (47) gives 


X c x c 
s s 


1_ 

B 


77 \Vi - Vi - 1, 


o ^ < 1 


(50) 
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0< £<i 


(51) 



1_ 

B 


Va - Vi + i> 


It is worth noting that the normalized coordinates X g , Y g of the interface yield a single 
interface curve that is independent of any parameter. The dimensionless coordinates 
x g /b and y g /b depend only on the physical parameter B. 

A relation is needed for heat flow through the frozen layer. Because the walls that 
form the corner are semi-infinite, the total heat flow through each wall depends on the 
total dista n ce along the wall from the corner. Since the frozen layer configuration is 
symmetrical about the corner angle bisector, only one wall need be considered. The 
heat flow through the vertical wall can be found by using equations (36) and (37). When 
c = 1 and d = 1, the denominator of equation (36) gives 


lim 
c— 1 
d—1 


r d{ — 

A V (| 2 -m + C)(d- £) 


lim 

dll + 1 )( c + *) 


K 



(c - l)(d - 1) 
(c + l)(d + 1)_ 


= K(0) = - 
2 


Then, the normalized heat flow through the vertical wall from the corner to an arbitrary 
location is given by 


Q 

4k(t f - t w ) 



0 < £ < 1 


(52) 


To find the value of Q/4k(t f - t w ) for a quadrant of a rectangular duct, equation (52) is 
applied to each of the two sides that represent the duct walls, giving 


Q 

4k(t f - t w ) 



^ V-*a 



(53) 


where £ and £, are found from equation (51) applied along each of the sides 
^a b 
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b b \ • Vi- v^/ 

B \ » Vi-V^/ 

Evaluation of Heat Flow From Interface Length 

As an alternative to obtaining the heat flow by the analytical expressions given pre- 
viously, a simple graphical method utilizing the interface curves can be used. If s is 
the entire length of the interface, then 

Q = h z s(tj - t f ) 

or, in dimensionless form, 

Q - h * b s V^ IMb (54) 

4k(t f -t w ) k 4bt f -t w A\b) 

For one quadrant of the duct, the interface length is s/4 and can be measured graphi- 
cally. This measurement is then divided by the length of the vertical side and multiplied 
by the parameter B to give the dimensionless heat flow. 

Some of the analytical expressions such as equation (53) are difficult to evaluate ac- 
curately in some instances. For aspect ratios larger than about 4, the £ a and £ b can 
become close to 1 so that the natural log term is very sensitive to small inaccuracies in 
4 a and In this instance, equation (54) is especially useful. 

For a very thin layer, the length of the interface is closely approximated by 
4a + 4b - 8y, where y is the frozen layer thickness when the layer is very thin relative 
to the duct dimensions . Substitute into equation (54) to obtain for thin layers 

9 -A /& — + 4 - 8 — \ B = /— + 1 - — \ B (55) 

4k(t f - t w ) 4 V b b j \b B ) 
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RESULTS AND DISCUSSION 


In the ANALYSIS were derived relations that can be used to evaluate the coordinates 
of the frozen interface and the heat flow through the frozen region as a function of the 
imposed physical conditions. This evaluation can be made for each aspect ratio. The 
desired quantities are interrelated by means of the mapping quantities c and d. The 
general analysis for c > 1, d > 1 was used first to evaluate frozen profiles by starting 
with large values of c and then decreasing c toward 1. The large values of c and d 
correspond to thick frozen regions. For a fixed aspect ratio, both the c and d de- 
crease as the frozen region becomes thinner. When c is very close to 1, some of the 
elliptic integrals become difficult to evaluate accurately. The limiting solution for 
c = 1 and d > 1 was then used by starting with the smallest value of d used in the 
general solution. Profiles of the interface were evaluated for several smaller d values 
until d was very close to 1. As d approaches 1, it is convenient to use the limiting 
case for c = 1, d = 1 applicable for thin frozen layers. 

The contours of the frozen region for aspect ratios a/b = 1, 2, 3, 4, and 5 are given 
in figures 9(a) to (e). The curves in each plot are for various values of the controllable 
physical cooling parameter (k/h^b)^ - t w )/(t^ - t^)]. The dimensionless heat flow as a 
function of this physical parameter is given in figure 10 on both rectangular and semi- 
logarithmic coordinates, the latter to better reveal the details for thick frozen regions. 

As a typical set of frozen configurations, examine figure 9(d) for a/b = 4. For thin 
layers, corresponding to a small value of the cooling parameter, the layer is of constant 
thickness except close to the corner. As the cooling parameter is increased within the 
range of (k/hjb)[(tj - t w )/(t^ - t^)J up to about 0.4, the frozen thickness increases fairly 
uniformly around the duct. Then, as the cooling parameter is further increased by a 
very small amount in the neighborhood of 0.5, the thickness along the short side in- 
creases substantially while that on the long side remains almost constant. For thick 
frozen regions, the interface approaches a circular shape. 

One of the interesting aspects of the profiles in figure 9 is that for each aspect ratio 
there is a maximum value of (k/h^b)^ - t w )/(t^ - tj)j. In the case of a/b = 4, for ex- 
ample, the maximum value is about 0.5. For a larger cooling parameter than this maxi- 
mum, no equilibrium profiles are possible, and the channel would freeze completely. 

The fact that there is a maximum (k/h^b)^ - t w )/(t^ - t^)J for each a/b leads to 
the most interesting feature of figure 9. For each value of the cooling parameter below 
the maximum, there are mathematically two possible frozen regions . One set, shown in 
figure 9 as solid curves, corresponds to thin and moderately thick frozen regions, while 
the dashed set is for very thick regions. Since the heat flow, as given by equation (54), 
is proportional to the interface length, the heat flow will also be a double- valued function 
of the physical cooling parameter. This function is shown in figure 10. As the channel 
aspect ratio increases, the heat flow (and interface profiles) are very sensitive to the 
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(k/hjb)[(tj - t w )/(tj - tj)J near its maximum value. As will be shown later in this sec- 
tion, the dashed set represents unstable contours that will not exist in actual physical 
situations. 

It is convenient to examine the simpler one- dimensional situation of freezing inside 
a cylindrical tube to discuss more easily the stability of the dashed and solid profiles . 

As shown in figure 11, if the inside radius of the frozen region is r^ the energy con- 
vened to the interface is Q conv = h^(t^ - tj) 27 rr.. This expression is equated to the con' 
duction through the frozen region to give 


t f - t 

" V 27 ^ = 

In — 

r i 

Rearranging yields 

k t r t w . r i lj r ° 

Vo h - *£ r o r i 


(56) 


(57) 


Equation (57) is plotted in figure 11, and the double- valued nature of r./r Q for various 
values of the physical parameter (k/h^r 0 )[(tj - t^)/^ - t^)] is evident. The maximum 
value of the cooling parameter is found by differentiating equation (57) and setting it 
equal to zero 


d 



= 0 


-1 + In — 
r i 


This operation gives r^/r Q = 1/e, which is then substituted into equation (57) to give 
the maximum value of the cooling parameter as 


k *f " *w 

h.r t. - t f 

l o l f 


max 


1 

e 


For a cooling parameter larger than 1/e, there is no condition where the convective heat- 
ing from the warm liquid can be as large as the conduction away from the frozen inter- 
face; for this condition, the tube will freeze solid. 
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The curve in figure 11 represents a state of equality between the convective heat 
supplied to the frozen layer and the heat conducted through the frozen layer. For every 
value of (k/h^r Q )[(t f - t w )/(t^ - t f )] < 1/e, two frozen-layer thicknesses are possible. 

The question then arises as to whether both layer thicknesses can exist as states of 
stable equilibrium. The answer is found by examining changes that occur when each 
equality state is perturbed slightly; that is, the equality radius ratio r./r Q of the frozen 
layer is increased or decreased slightly and then examined to see if the radius returns to 
its original value. 

The stability analysis is begun by noting that Q conv changes linearly with rj and 

the conduction through the layer Q cond varies as [ln^/r^] -1 . The difference Q conv 
- Q cond is from equation (56) 


Qconv - Qcond “ 2 "o(tWz - V - 2 ” k tf ’ ‘ W 

Vo/ 



(58) 


At equilibrium, Q conv - Q cond = 0 at every r^ Now, examine the small change 

d(Q - Q j) that occurs when a small disturbance causes a change dr. in the loca- 
conv cona ^ 

tion of the liquid - frozen-layer interface. To do this, write 



Differentiating equation (58) with respect to r./r Q results in 



Since the interface is being perturbed from equilibrium, equation (57) is used to simplify 
the expression in brackets to obtain 
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d 


(59) 


lllllllllllllllllllllllllllllll 



(Q, 


conv 


- Q cond> - 2vr o - V 


1 +• 



When r^/r o > 1/e, the derivative is negative, and, for r./r Q < 1/e, the derivative is 
positive. 

Consider first the region 1 > r^/r Q > 1/e, which is along the solid portion of the 
curve in figure 11 and for which [d/d( r^/r Q )j (Q c Qnv - Q cond ) is negative. If r./r Q is 
increased by d(r./r ), Q conv increases because of the area increase at the interface, 
and Q cond also increases because of the decreased frozen- layer thickness. However, 
d(Q CO nv " Q CO nd^ * s ne S a ^ ive in this region, which means that as r./r Q increased, 

Q CO nv increased less than Q conc j increased. That is, the cooling by conduction in- 
creased more than the heat supplied by convection. Hence, refreezing occurs, and the 
original radius is restored. With similar reasoning, when r./r Q is decreased by 
d( r i/ r ), ^conv decreases because of the area decrease, and Q cond decreases because 
of the increased frozen-layer thickness. However, because d(Q conv - Q CO nd^ is P osi " 

tive, Q„ dominates Q ,, and the added increment d(r./r ) of frozen layer is 

melted to restore the original radius. This establishes that the solid portion of the 
curve in figure 11 gives states that are in stable equilibrium. 

Now, consider the region 0 < r^/r < 1/e, which is along the dashed portion of the 
curve in figure 11 and where |d/ d( r .. / r Q )J (Q conv - Q con( j) is positive. When the perturba- 
tion d(r./r Q ) is positive, the d(Q conv - Q cond ) is also positive; consequently, Q conv 
increases more than Q cond - The frozen layer is unstable and will melt until the stable 
equilibrium radius on the solid portion of the curve is reached (path a in fig. 11). 

When d(r i /r Q ) is negative, d(Q conv - Q cond ) is negative; therefore, the increase in 
Q con d dominates, and freezing occurs until all the liquid is frozen (path b in fig. 11). 
Thus, the dashed portion of the curve in figure 11 represents an unstable equilibrium 
region that will not actually exist in reality. 

In summary, for the circular tube case for the situation where a flow of warm liquid 
is maintained, the only stable configurations are when 


k 


h.r 
l o 



and 



r o e 


For(k/h z r 0 )[(t f -t w )/(t r t f )] > 1/e, no equilibrium states exist, and the tube will freeze 
completely. The tube can also freeze completely for (k/h^r Q )[(t£- t w )/(t^ - t^)J < 1/e if 
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by some means the r^/r Q falls below the dashed curve in figure 11. 

This same type of reasoning applies for the rectangular ducts. The maximum val- 
ues of (k/h^b) [(tj - t w )/(t^ - tj)], for which frozen regions can be maintained with flowing 
liquid, can be found from figure 10. The values ranged from about 0.4 to 0. 5 as the as- 
pect ratio is increased from 1 to 5. The solid profiles in figure 9 are the stable pro- 
files. The dashed profiles are admissible in the sense of a mathematical equality of heat 
flows but are unstable physically. 


CONCLUDING REMARKS 

A conformal mapping technique was utilized to determine two-dimensional frozen- 
layer configurations for solidification in a cooled rectangular channel carrying flowing 
warm liquid. The frozen region was mapped onto a potential plane W = -T + ii}/, where 
the negative of the dimensionless temperature is the potential function for heat flow. 

The shape of the frozen region is found from the integral relating the physical plane to 
the potential plane 

Z = / — L_ dW 

J C(W) 

where £ = -(9T/3X) + i(3T/3Y). The functional relation between £ and W, necessary to 
carry out the integral, was found by conformal mapping. 

^ The configuration of the frozen region for each duct aspect ratio depended on a sin- 
gle physical cooling parameter (k/h^b) [(t^ - t w )/(t^ - t^)J. For each aspect ratio, there 
was a maximum value of the cooling parameter above which the duct would freeze com- 
pletely. The interesting feature found was that below this maximum, two possible mathe- 
matical solutions exist for each value of the cooling, one was a thin frozen region and the 
other a thick region. Stability considerations revealed, however, that the thick region 
configuration was unstable. As a consequence, only thin frozen region configurations 
are physically possible and, then, only when the value of (k/h^b) j(tj - t w )/(t ^ - t^)] is 
below a certain magnitude for each aspect ratio. 
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(a) Potential plane, W = -T + 'up. 




(c) Intermediate w plane, o) = \n\c\ + i0. 
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td) Intermediate u plane, u - £ + \n. 

Figure 3. - Transformation planes for general case: c > 1, d > 1. 
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(a) Layer is one-dimensional in central region of long sides, (a) Dimensional coordinates in z plane, 

c - 1, d > 1. 




(b) Layer is one-dimensional in central region of long and 
short sides, c - 1, d - 1. 

Figure 4. - Limiting cases where portions of frozen region 
are one-dimensional. 
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(b) Normalized coordinates in Z plane. 

Figure 5. - Limiting case of frozen layer one 
dimensional along long side of channel. 



Aip 

3 


2 

4 


1 

(-1, 0) 


7 -T 

3 

T = 

■ 1 

6 


T = 0 


(a) Potential plane, W = -T + i ifi. 



(b) Temperature derivative plane, 
£ = -(dT/dX) + KdT/dY). 
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(c) First intermediate plane, (d) Second intermediate plane, 

oj = ln|£| + i0. u = £ + fy. 

Figure 6. - Transformation planes for infinite aspect ratio duct. 
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(a) Dimensional coordinates. (b) Normalized coordinates. 

Figure 7. - Limiting case of very thin frozen layer in channel corners. 



(a) Potential plane, W = -T + i ip. (b) Temperature derivative plane, 

£ 2 -(dT/dX) + i(dT/dY). 



(c) Intermediate u) plane, (d) Intermediate u plane, 

oj= ln|£| + ie. u = £ + bj. 


Figure 8. - Transformation planes for corner solution. 
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Thin frozen region, stable 



fa) Rectangular coordinates. 
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Cooling parameter, (k/hjb)[(t f -t w )/{tj-t f t] 

(b) Semilogarithmic coordinates. 

Figure 10. - Relation between total heat transfer and parameter con- 
taining imposed conditions. 






Figure 11. - State of equality between heat convection 
from liquid and heat conduction in frozen layer 
inside a circular tube. 
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